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1.0 SUMMARY 

Three types of navigation onorbit numerical integrators were evaluated, and the 

following results were obtained: 

a. Power integrators with no delta-V incorporation, just coasting; i.e., using 

Taylor series expansion integrators 

(1) Super G is slightly better than average G for step sizes of 2 and 4 
seconds. (A 1200-meter error for delta-T = 4 seconds after 10 revolu- 
tions (revs)). Super G is marginal for delta-T = 15 seconds. Neither 
are adequate for delta-T > 30 seconds. 

(2) Spiffy G shows no improvement over super G. 

(3) Super G4 (third order) has slight improvement over super G or spiffy G 
at steps of 2 and 4 seconds, but is inferior to them for delta-T >15 
seconds . 

b. Coasting integrators using the Cowell method of special perturbations 

(1) With the exception of Runge-Kutta third order, all third-order (RK or 
Nystrom) integrators performed rather poorly for delta-T >15 seconds. 
The RK3 is a remarkable exception, and it competes favorably with 
fourth-order integrators with delta-T up to 60 seconds. (The RK3 error 
is less than 1000 meters for delta-T = 60 seconds and 10 revs). 

(2) The fourth-order Nystrom integrators performed as well, or slightly 
better than the RK's but they are a little slower to exeerte. The 
Nystrom 4 with Lear's coefficients ((ref. 1) - NLXD4/4) performed 
better than all other Nystrom integrators. 

(3) All fourth order integrators at delta-T = 2 seconds had a 0.1 -meter or 
less error when compared to the KS (ref. 2) reference integrator. 

(4) In general, RK4 integrators perform adequately for up to AT = 60 sec- 
onds and 10 revs with errors less than 1200 meters except for RKL42. 
Degradation occurs rapidly beyond that AT with the exception of 
RKL41 (ref. 1), which is adequate for up to 120-second time steps. 

c. Coasting integrator using the Pines variation of parameter perturbation 

method . 

(1) The Pines formulation with RKG4 (the standard predictor for onorbit nav 
igation, ref. 3) performs excellent for up to 5-minute (300 seconds) 
time steps (error less than 200 meters for delta-T = 5 minutes and 10 
revs). The 10-minute time steps may be considered. 

(2) The Pines method used more core and executes slower than the Cowell 
method for a single step. However, for certain applications where 
longer time steps are permitted, this method is more time efficient. 
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2.0 INTRODUCTION 


This document presents results for using the three onorbit navigation integrators 
in the onboard software: (a) average G for user parameter propagator (UPP), (b) 

super G for the onorbit navigation state propagation function, (c) Pines/RKG for 
the onorbit state prediction function, and (d) other potentially useful Runge- 
Kutta and Nystrom integrators for onorbit navigation where an analysis task was 
performed with typical Shuttle orbits (i.e., 100- to 300-n. mi. altitude and 
small eccentricities). 

The acceleration function for a simulated force model (app. E) included a cen- 
tral force field, J2 gravity terms, and a drag perturbation (ref. 4). These 
were programed into a Hewlett Paokard HP9825 desktop calculator (12-digit ma- 
chine; no double precision). Gravity up to 4x4 (fourth degree-fourth order) was 
also investigated by CSDL, and its results are included here for completeness. 

All oases were run for approximately 10 revs or until position error (RSS) was 
greater than 100 kilometers. 

Integrator step sizes considered were 2, 4, 15, 30, 60, 150, 300, and 600 
seconds . 

The following integrators were considered in this analysis for total accelera- 
tion integration (ref. 1): 

a. Series expansion integrators for powered flight: 

(1) Average G - second order 

(2) Super G - second order (app. A) 

(3) Spiff y G - second order (app. A) 

(4) Super G4 - third order (app. A) 

b. The RK/NYSTROM integrators using the Cowell method for coasting flight: 

(5) RK3 - standard Runge-Kutta third order 

(6) RKL3 - Lear’s coefficient for RK3 

(7) NLXD4/3 - Nystrom-Lear coefficient for NXD4/3 

(8) NXD4/3 - Nystrom third order 

(9) RKG4 - Runge-Kutta-Gill fourth order 

(10) RKL41 - Runge-Kutta-Lear fourth order 

(11) RKL42 - Runge-Kut* i-Lear fourth order 

(12) RK4 - standard Runge-Kutta fourth order 
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(13) NLXD4/4 - Nystrom-Lear coefficient fourth order 

(14) NXD4/4 - Nystrom fourth order 

o. The following variation of parameters for special perturbation was examined 
for the ooasting flight: 

(15) Pines /RKG4 (ref. 3) 

Three integrator initial conditions (IC) were used for this analysis and are 
listed in table I. The positions and velocities (in kilometers and km/ sec) 
listed constitute the state vector at time = zero second. After approximately 
10 revolutions (t = 54000 seconds), the state vector propagated by the reference 
integrator - a KS (Kustannheimo-Stieffel) formulation (ref. 2) with Runge-Kutta 
45 integrator (table I) was compared with the tested integrator. The position 
difference was RSS (root sum squared) to determine an integrator error for the 
various integrator steps attempted. The position errors are listed in tables II 
through V. 

Integrators (2) and (15) (super G and Pines/RKG) were also tested by Charles 
Stark Draper Laboratory (C5DL) personnel in the HAL code environment (ref. 5). 
The CSDL results basically duplicated the results of tables II to V. In addi- 
tion, reference 5 provides actual execution time for API 01 (Shuttle onboard 
computer) in extended precision for 1, 5, and 10 revolutions (revs) of propaga- 
tion for various step sizes. Some of the data in reference 5 are duplicated in 
this report. 


3.0 ANALYSIS RESOLTS 


3.1 CASE 1 DATA 

For case 1, all 15 integrators were ini i*Uzed with IGfl and allowed to run for 
approximately 10 revs (54 000 seconds) with state vectors printed at 10-minute 
increments. Only J2 gravity perturbation was included in the functional evalua- 
tion call to the acceleration function. (Drag was used in case 4 only.) ICfl, 
which is a 146-n. mi. circular orbit inclined at 30° with the equator, was used 
as the basic orbit to screen out integrator performance. The integrators that 
performed well in this environment were further evaluated in cases 2, 3, and 4. 

To quickly assess integrator performance, the energy and percent delta energy 
equations were programed. Although no energy data are presented in this re- 
port, it was found that with these parameters, coding errors were detected much 
sooner than by the normal differencing of state vectors along the reference 
trajectory. The initial orbit energy Kq was printed at the beginning of the 
run and subsequently, delta energy AC was printed. For conservative orbits; 
i.e., no drag or self-induced satellite accelerations like uncoupled RCS thrust, 
the delta energy must be zero along the trajectory: i.e., energy is conserved. 
(This is a necessary but insufficient condition that indicates to the user that 
orbit errors are probably not being introduced by the integration scheme 
selected.) It was found that with 1 digits of delta energy printed, integrator 
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induced errors in the orbit are detected ouch sooner (in 10 or 20 steps) than 
they would by differencing the reference state vector with the tested integrator 
state vector. Therefore, integrator errors could be detected within a few 
time steps rather than after a rev of data. 

The following formulas were used for energy and delta energy: 



A 




Ci 


P 

R 


V2 

Z 


A 


ACi = 


Ci - Co 

Co 


where 


A 

= orbit energy at time = t^ 


A 

= initial orbit energy at time = 0 


p = Earth gravitational constant = 398601.0 km^/sec 2 


A 

R = satellite position vector magnitude, km 


A 

V = satellite velocity vector magnitude, km/sec 
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A s potential function 


A 

X = J2 perturbation constant 


A 

X.j , Xg, = inertial components of satellite position vector, 


X-| along Earth's equator, X3 along Earth's North Pole 


Kj = J2 potential constant s 1.08265 x 10*3 


A 

R e = Earth's radius = 6371.22 km 


A 

A£i = delta energy from time = 0 


Table II shows that the average G integrator performs adequately for steps up to 
4 seconds and then degrades rapidly. Super G performs slightly better and is ef- 
fective up to 15-second step sizes before collapsing. Spiffy G (app. A) 
performs almost identically to super G in all cases. As mentioned in section 
1.0, a 1200 meter error (fig. 1) was noted for super G after 10 revs of propaga- 
tion for 4-second time steps. 

Super G4, developed in appendix A (ref. 6) performed quite well up to 4-second 
step sizes but degraded quickly for higher steps and was not evaluated further. 

The fourth-order Runge-Kuttas used a common fourth-order RK algorithm and only 
the coefficients were changed for each integrator. Functional flow charts for 
the third- and fourth-order Runge-Kutta and Nystrom integrators are given in 
appendixes B and C. All coefficients were obtained from reference 1 and are 
listed in appendix 0. 

The functional evaluation subroutine was obtained from reference 7 and is shown 
in the appendix E flow chart. This flow chart was used by integrators to deter- 
mine the acceleration vector. Note that only central force field, J2, and drag 
is used. 

All fourth-order Runge Kutta integrators (RKG4, RKL41, RKL42 , and RK4) performed 
quite well with delta steps up to 60 seconds. The RKL41, a Runge-Kutta 
integrator with Lear's first set of coefficients performed quite well with delta 
steps up to 150 seconds. 
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TWo third-order Runge-Kutta integrators were tested (a standard RK3 and RKL3 
that used the Lear coefficients). The RK3 performed exceptionally well for a 
third-order integrator and, in fact, performed almost as well or better than the 
fourth-order Runge-Kutta integrator for delta-T < 60 seconds. The RKL3 did not 
perform well. 

Two fourth-order Nystrom integrators were tested: the NLXD^ and the NXD4. TT~ 

NLXD4, which used the Lear coefficients, performed slightly better than the st • 
dard Nystrom NXD4. In general, the Nystrom integrators performed at least as 
well as the RK's. However, the algorithm took slightly more time to execute be- 
cause of the extra calculations required in the Nystrom algorithm. 

Two third-order Nystrom integrators were tested, and as shown in table II, they 
performed rather poorly and were quickly discarded (NXD3 and NLXD3). 

The last integrator tested was the fourth-order Runge-Kutta -Gill (RKG4) using 
the Pines variation of parameters formulation technique. The code was obtained 
directly from the onorbit navigation FSSR (ref. 3). As shown in table II, this 
formulation and integrator combination performed exceptionally well up to 10- 
minute steps. However, the Pines code is more complex than Cowell's formula- 
tion, and as shown in table VI, it does take about five times longer to execute. 

Table VI lists the relative time it took for the various algorithms to perform 
?n the HP9825 environment (a fairly accurate but somewhat slow desktop calcula- 
tor when compared with a typical powerful and much faster machine such as the 
UNIVAC 1108). It should be noted that the algorithms coded were general purpose 
and inefficient, specifically in terms of execution time because if a coeffi- 
cient of zero were encountered, the algorithm operation was still performed. In 
any case, this table clearly shows that Pines/RKG is not only the most accurate 
method but also the slowest method tested for a single step. More or less the 
same relative execution timing data were observed in reference 5 and are 
reproduced elsewhere in this report. 


3-2 CASE 2 DATA 

For case 2, the initial condition IC#2 represented a Skylab reboot rendezvous 
orbit after the terminal phase finalization (TPF) maneuver. To conserve com- 
puter time, some of the small delta step runs were deleted. It was assumed that 
the trend to higher accuracy for smaller time steps had been established in the 
case 1 results. Only the average G, super G, spiffy G, RKG4, RKL41, RK3, and 
the NLXD4 integrators of case 1 were tested for this case. 

Results for case 2 were basically the same as for case 1 and are tabulated in 
table III. 


3.3 CASE 3 DATA 

Case 3 results are listed in table IV. The IC#3 represented a Skylab insertion 
orbit of 100 n. mi. circular. Because drag is significant at this altitude, the 


6 



79FM25 


integrators of oase 2 were tested with the drag perturbation of reference 4 (a 
simplified atmospheric model) in the acceleration function. 

The accuracy results for case 3 were very similar to those of cases 1 and 2. 

This seems to indicate that drag, if modeled properly in the acceleration func- 
tional evaluation used by the integrator, will present no difficulty to the Inte- 
gration scheme used. However, it was noted immediately that the orbit energy 
and delta energy computations were being affected by the slight acceleration 
produced by the drag. Therefore, if drag or some other nonconservative force 
are included in the acceleration model, the value of the energy check in 
evaluating integrator performance is lost. 


3.4 CASE 4 DATA 

Case 4 results are given in table V. In essence, this case is ident.o- 1 to case 
3 except that no drag was included in the evaluation of the acceleration. _->me 
runs with super G4 and Pines/RK6 were added to get more data for this integra- 
tor. As with previous cases, the accuracy results were very similar to the 
other cases. The integrator error differences between cases 3 and 4 (i.e., drag 
versus no drag) were very minimal. However, the actual position differences 
after 10 revs were about 30 kilometers (for constant Shuttle surface area) due 
to the drag. 


3.5 CSDL RESULTS 

To obtain some performance data for OPS-2 state propagators in HAL code running 
in AP101 extended precision, CSDL used IC#1 and IC#2 and propagated them for 10 
revs for both the super G and Pines/ RKG integrators at various integrator step 
sizes in their API 01 simulation. 

Since the onboard flight code was to be used for this analysis, it was decided 
that an existing ENCKE-Nystrom formulation would be used as a real world refer- 
ence test case. This formulation consists of a full double precision with a 
fourth-order /degree gravity model and an integration step size of 30 seconds. 
(This was later changed to an 8 degree/order gravity model to determine the dif- 
ferences between a 4/4 and 8/8 gravity). Results are given in reference 
5 and they can be summarized as follows: 

a. For J2 only (gravity s 2 order, zero degree), which is what the HP9825 pro- 

gram was formulated to, the differences between the super G and RKG/Pines 
were about 1000 meters - basically the same as tables II and V. Reference 5 
results can be summarized in table VII. (They can be computed by 
differencing cases: super G6 and PRKG5; super G6' and PRKG5'.) This same 

integrator difference is basically observed for cases usi.g higher fidelity 
gravity models: i.e., super G1 and PRKG1; super Go and PRKG3; super G5 and 

PRKG4; and the similar primed cases. 

b. Table VII shows that the Pines/RKG, using a 3-minute step size (predictor 
type operation), performs almost the same as for the 1 -minute step size. 
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4.0 INCLUSIONS AND RECOMMENDATIONS 

a. The super G integrator is a very simple and effective integrator for 2- and 
4-seoond time steps. Since IMU delta-V data oan be easily incorporated in 
the integration soheme , its use as the standard onorbit navigation 
propagator for the maintenance of the current state has been implemented in 
the onboard navigation software. 

b. The Pines variation-of-parameters formulation method with a Runge-Kutta-Gill 
(RKG) fourth-order integrator method pro-: ices excellent results up to 300- 
second time steps. Onorbit prediction w~wh this method (3- to 5-minute time 
steps) has been implemented in the onboard onorbit navigation scheme. 

c. The Runge-Kutta third order (ref. 1), using Cowell's method, is an excellent 
general purpose orbit determination integrator for time steps up to a 60- 
second duration. 
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TABU I.- INTEGRATOR INITIAL CONDITIONS 


Param 

Unit 

Case 1 

Case 2 

Case 3 

Case 4 

t c = 0.0 

*o 

km 

6 649.02 

-3 972.220046 

-4 192.451762 

-4 192.451762 

*o 

km 

0 

-1 892.121621 

4 777.342541 

4 777.342541 

Zo 

km 

0 

-5 000.635973 

1 636.711626 

1 636.711626 

*o 

km 

0 

4.393527408 

-4.85227183 

-4.35227183 

*o 

km/sec 

6.705343087 

-6.262423238 

-2.327477093 

-2.327477093 

K 

km/sec 

3.871331637 

-1.112883779 

-5.64056083 

-5.64056083 

»A 

km 


284.08 

185.05 

185.05 

Hp 

km 


269.85 

181.38 

181.38 

i 

deg 

30° 

49.86° 

49.93° 

49.93° 

Final position 6 tp = 54 000 sec 



W/0 drag 

W/0 drag 

W drag 

W/0 drag 

*F 

km 

6 507.6213 

-4 068.7385 

-4 955.392 

-4 968.7873 

*p 

km 

1 027.5009 

-1 666.3453 

- 543.481 

- 519.5312 

Z F 

km 

895.9049 

-5 003.6163 

-4 263.2106 

-4 211.5175 
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TABLE II.- CASE 1 


Integrator 

2 sec 

4 sec 

15 sec 

AVE G 

759 

3035 

42 687 

SUPER G 

332.2 

1 161.1 

3 360.8 

SPIFFY G 

332.5 

1 161.2 

3 360.8 

SUPER G4 

22.9 

37.2 

8 857.3 

RKG4 


.1 

1.1 

RKL41 


.1 

.1 

RKL42 

.1 

.1 

10.5 

RK4 

.1 

. 1 

1.8 

RK3 

.1 

.3 

9.8 

RKL3 

21.3 

169.6 

8 939.8 

NLXD4 


.1 

.1 

NXD4 

.1 

.1 

2.0 

NLXD3 

150.7 

918.6 

17 747.1 

NXD3 

390.7 

2 406.8 

48 173-5 


PINES-RKG4 


.1 
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170 

682 


681 

716 






57 

526 


228 

962 






57 

526 


228 

962 






70 

776 


563 

902 







16 

.3 


227 

.1 

4 

558.5 

92 

327.1 



.8 



.7 

1358.3 

50 

861 


244, 

.5 

6 

326 


515 

726 




43 

.5 

1 

159 

.7 

100 

114 




66, 

.9 


164 


84 

411.7 



71 

545 

.3 









1 , 

,0 


31 

.7 

3 

121.6 

103 

271 


34, 

.8 


672 

.1 

39 

134.4 



82 

921 , 

.5 









233 569 


.1 


.1 


1.70 


33.4 


683.9 
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TABLE III.- CASE 2 INTEGRATOR POSITION ERROR AFTER 10 REVS, METERS 

No drag; i = 49.8°" 

Skylab TPF 

H a = 153.4 n. mi. 

Hp = 145.7 n . mi. 






Delta-T 



Integrator 

2 sec 

4 sec 

15 sec 

30 sec 60 sec 

150 sec 

300 sec 

AVE G 

752.3 

3 009.4 

42 312 

169 187 675 543 



SUPER G 

328.7 

1 149.7 

3 434.9 

55 915 783 477 



SPIFFY G 



3 434.9 

55 915 783 477 



RKG4 


.1 

1.0 

16.0 223.3 

4 500.5 

89 065.1 

RKL41 



.1 

.8 2.2 

1 228.4 

47 179 

RK3 



.6 

10.8 751.0 

91 251 


NLXD4 



.1 

1.0 30.9 

3 041.6 

100 600 
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TABLE IV.- CASE 3 INTEGRATOR POSITION ERROR AFTER 10 REVS, METERS 

With drag; i * 49,93°"| 

Skylab insertion J 

100-n. mi. circular I 


Integrator 

2 see 

4 sec 

15 sec 

30 sec 

60 sec 

150 sec 

300 sec 

AVE G 

782.9 

3 133.4 

44 064 

176 119 




SUPER G 

349.0 

1 215.4 

3 044.3 

64 646 

878 152 



SPIFFY G 

349.5 

1 216.3 

3 044.2 

64 659 

878 232 



RKG4 


.1 

1.1 

15.4 

209.0 

3 526.5 

126 550 

RKL41 


. 1 

. 1 

.8 

.7 

1 487 

54 573 

RK3 


.1 

3.9 

14.8 

"37.6 

103 727 


NLXD4 


.1 

.1 

1.2 

35.6 

3 501.4 

115 984 



7W^ 




table 

V.- CASE 4 

INTEGRATOR position error 

[No drag; i = 49.93° 

I Sky lab insertion 
[100-n. mi. circular 

AFTER 10 REVS 

t METERS 



Integrator 

2 sec 

4 S»*C 

15 sec 

jO sec 

60 sec 

150 sec 

300 sec 

600 sec 

1200 sec 

AVE G 

77*».9 

3 100.1 

43 591 

174 299 

695 897 





SUPER G 

347.3 

1 207.9 

3 016.8 

64 267.8 

871 275 





SPIFF* G 

347.2 

1 208.0 

3 016.8 

64 267.8 

871 275 





SUPER G4 

46.1 

269.1 








RKG4 


.1 

1.1 

17.1 

237.7 

4 427.6 

116 140 



RKL41 


. 1 

.2 

.4 

7.3 

1 775.3 

59 263 



RK3 


.1 

3.7 

14.3 

744.7 

103 285 




NLXD4 


.1 

0 

1.0 

35.2 

3 478.7 

115 263 



PINES-RKG4 



.2 

.1 

.1 

5.4 

122.7 

1 247.7 

42 749 
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TABLE VI.- 

INTEGRATOR RELATIVE TIMING DATA FOR 
bP9825 EXECUTION TIME 

Integrator 

(a) 

Order 

Approximate time, 
step/sec 
'b) 

Average G 

2 

0.293 

Super G 

2 

.320 

Spiffy G 

2 

.327 

Super 04 

3 

.643° 

RK3 

3 

.493 

RKL3 

3 

.493 

NLXD4/3 

3 

.527 

NXD4/3 

4 

.52? 

RKG4 

4 

.693 

RKL41 

4 

.693 

RKL42 

4 

.693 

RK4 

4 

.693 

NLXD4/4 

4 

.720 

NXD4 

4 

.720 

Pines 

4 

3.375 


a Nonoptimum code could be reduced significantly 
for some, especially super G4. 
including F evaluation J2; no drag. 
°Inefficient coding by programer. 
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TABLE VII.- CSDL RESULTS 


Case 

Time, 

Gravity 

Model 

IC 


Total position 

error, ft 

no. 

step/sec 

deg 

order 

no. 

Integrator 

1 rev 

5 revs 

10 revs 

1 

4 

2 

0 

i 

Super G 

2588 

13 650 

35 358 

3 

4 

2 

0 

2 

Super G 

2593 

6 064 

20 807 

2 

4 

4 


1 

Super G 

479 

2 182 

3 814 

4 

4 

4 

4 

2 

Super G 

475 

2 156 

3 771 

2k 

8 

4 

4 

1 

Super G 

1875 

7 623 

10 836 

4A 

8 

4 

4 

2 

Super G 

1856 

7 545 

10 761 

5 

60 

2 

0 

i 

Pines-RKG 

2134 

11 634 

32 000 

8 

60 

2 

0 

2 

Pines-RKG 

2137 

4 096 

17 714 

6 

60 

4 

4 

1 

Pines-RKG 

5 

160 

556 

9 

60 

4 

4 

2 

Pines-RKG 

9.5 

1 61 

609 

7 

180 

4 

4 

1 

Pines-RKG 

19 

185 

540 

10 

180 

4 

4 

2 

Pines-RKG 

33 

415 

1 499 

6A 

60 

4 

2 

1 

Pines-RKG 

473 

1 602 

10 204 

9k 

60 

4 

2 

2 

Pines-RKG 

2655 

8 248 

21 888 

6B 

60 

2 

2 

1 

Pines-RKG 

390 

2 751 

13 322 

9B 

60 

2 

2 

2 

Pines-RKG 

2126 

5 544 

16 543 

2’ 

4 

4 

2 

1 

Super G 

156 

3 469 

13 383 

4' 

4 

4 

2 

2 

Super G 

3117 

10 193 

25 060 


Note: Truth model for above data used an Encke-Nystrom formulation in full 

double precision with a fourth-order and degree gravity model and an 
integration step size of 30 seconds. 
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TABLE VII.- CSDL RESULTS (Concluded) 


Case 

Time, 

Gravity 

Model 

IC 


Total position error, ft 

no. 

step/sec 

deg 

order 

no. 

Integrator 

1 rev 

5 revs 

10 revs 

2" 

4 

2 

2 

1 

Super G 

413 

4 698 

16 424 

4" 

4 

2 

2 

2 

Super G 

2592 

7 559 

19 780 

2A* 

8 

4 

2 

1 

Super G 

1430 

8 844 

20 402 

4A' 

8 

4 

2 

2 

Super G 

4494 

15 575 

32 043 


Note: Truth model for above data used an Encke-Nystrom formulation in full 

double precision with a fourth-order and degree gravity model and an 
integration step size of 30 seconds. 
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APPENDIX A 

SUPER G AND SUPER G4 EQUATIONS 


A-1 
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Super G and super G4 algorithms: 
a* Compute Fq ~ F(Tqi Rq> ^o^ 

b. Compute Ri , V-j 

Rl = R 0 + V 0 DT + F q DT 2 /2 
V 1 = V 0 + F 0 DT 

Ti = T Q + DT 

c. Evaluate Fi = F ( T i , Ri , 

d. Update R -j , V-j 

R 1 = Ro + v o DT + + A3 • DT 5 /6 

If super G *► Go to step (e) 

v ! = V 0 + F 0 DT + A3 * DT 2 /2 

Where 

A3 = (Ft - F 0 )/DT 

e. If (Spiffy G or Super G) Relabel R 0 R-| 

v o - V 1 

T 0 - T 1 

And go to (a) 

f. Evaluate F-j = F (T^, R^, V^) 

g. Compute R 2 , V 2 

R 2 = Ri + VtDT + F iDT 2 /2 + A3 * DT3/6 

DT 2 

V 2 = V-| + F tDT + A3 • — 

Where 



T 2 = T q + 2DT 


A-3 
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h. Evaluate F 2 = F ( Tg » R 2 > V 2 ) 

i. Update R 2 . ^2 

R 2 = Ri + V^DT + FiDT 2 /2 + A3 • DT 3 /6 
V 2 = V. + F^DT + A3 • DT 2 /2 
Where 

A3 = (F 2 - F^/DT 

j. Update F 2 = F (T 2 , R 2 , V 2 ) 

k. Update R 2 . v 2 

R 2 = Ri + V^DT + FiDT 2 /2 + A3 - DT 3 /6 + A4 • DT 4 /24 
V 2 = V! + F tDT + A3 * DT 2 /2 + A4 • DT3/6 
Where 

(F 2 - F o> 

A3 5 

2DT 

F 2 - 2Ft + F 0 

A4 

DT 2 

l. Evaluate F 2 = F (T 2 , R 2 , V 2 ) 

m. Relabel 

F o* F 1’ T o«* T 1 
R-| +■ R 2 , Vi «- v 2 , Fi «* F 2 , T-| *■ t 2 

GOTO (g) 
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Derivative a .iroximations : 

a. Two points (T Q , F 0 ) , (T ^ , F-|) Known 
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0(DT 2 ) 


The errors (DT) and (DT^) are due to the derivative approximations. 
Further errors may be introduced due to errors in F Q , F^, or F2* 
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APPENDIX B 

FUNCTIONAL " M CHART 
FOR RK3 AND R”U 
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Figure B1 . - RK3 functional flow chart. 



1 


CALL 

"ACC" 

( Vj • • 

• r 12» k 24» k 25> k 26^ 

Second functional evaluation 


J \ jr 

i 

i 

j - - 




For 

For 

i = 1 
i b 1 

to 3 ; 
to 6 ; 

k 2i = n(i + 9) ; next 
r( i 4 6) s ri + AT(ri6 

i 

• k u ♦ ri 7 • k2i ) ; next i 


i 

i 

1 




CALL 

"ACC” 

(ry . , 

i 

• r 12 » k 34» k 35 « k 36 j 

Third functional evaluation 


I 

I 


For i s 1 to 3 ; k 3 i s r(i ♦ 9) ; next i 

For i s 1 to 6 ; X* = r* + ATCr^g • k^ + * k 2 t ♦ r20 * k 3i) » next i 


< 3 - 


! 

! 

/ \ 

/ \ 

Mo / \ Yes i i 

n curh =\ PRINT: FINAL STATE I 

\ T f / ! AND STOP I 

\ / I ! 

\ / 

\ / 


Figure B1 . - Concluded. 



CD 

I 

cn 


C RK4 START 


LOAD COEFF: r 


t CURR * T i J J= 1 3 -* 25 r 


o-- 


RK4 LOOP 


I LOAD: 

*1 

* X 5; T T ^ j 

1 

at i 

I 

1 


Standard 

. RK4 coefficients shown. here 


1 

1 


r 15 8 1 

j r 13 » 

2 

r l*4 8 j 



1 



1 

; <" 16 * 

I 

O 

II 

t- 

u 


r l8 8 g 

• r i9 s 

0 

r 20 s 0 


r 21 8 1 

! 

1 


1 

1 

j r 22 8 

» - - 

6 

r 23 * rgi* 

= 3 

r? 5 8 6 


T C'JRR s t CURR + AT 
LOAD STATE VECTOR: 


r 1 .'X, r 2 = X 2 r 3 = X3 

r 4 s x 4 r 5 s x 5 r 6 s x 6 


POS 

VEL 


! input output I 

■ CALL "ACC" (r 1 . • . r$, k^, k^, k^) ! «• First functional evaluation 
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Figure B2.- RK4 functional flow chart. 





0 


CALL W ACC W (ry . . . k 24» k 25* *20 ‘ ♦ Second functional evaluation 


For i : 1 to 3 ; k 2 i * r(i + 9) ; next i 

?*** i : 1 to 6 ; r(l + 6) r r* + AT(ri 7 • + r^g • k 2 i) ; next i 


CALL n ACC w (Ry . . . r^ 2 , k^, k^5» k^g 


4 . Third functional evaluation 


For i= 1 to 3 ; = r(i + 9) ; next i 

For is 1 to 6 ; r(i + 6 ) = r* ♦ AT(rig ■ k- u ♦ r^ * k 2 i + r* 2 i • kj^ ) ; next i 


CALL "ACC" (ry . « , r -j 2 , kj^, ^ 45 » kjjg ) 


Fourth functional evaluation 


For i : 1 to 3 ; k ^ s r(i ♦ 9) ; next i 

For i s 1 to 6 ; X* r* ♦ AT(r22 m + r 23 * k 2i + r 24 * k 3i + r 25 # k 4i^ » next * 


/ \ 

/ \ 

No / \ Yes I 

/T CURR =\— PRINT: FINAL STATE VECTOR 

\ T f « / I AND STOP 

\ / I 

\ / 

\ / 


Figure B2.- Concluded. 
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APPENDIX C 

FUNCTIONAL FLOWCHART 
FOR NYSTROM FOURTH-ORDER INTEGRATOR 




0 

1 


>3 

M 

o 

w 

a 

»— » 

z 

o 

*13 

> 

o 

w 

CO 

r 

> 

Z 

* 

z 

o 

H 

t— t 

r 1 


PO 

8 5 

g £ 

c 5 

i» •** 

r O 

3 3 


(nYSTROM Hi STAKT^-- 


LOAD COEFF: rj !• 

t CURR 3 T i j * a 3 ♦ 35 I- 



■i NYSTROM 4 LOOP 


T CURR 3 T CURR + fiT 


LOAD: - X 6 ; T t ; T f ; 

1 

AT : 

j 

Standard NXD4/4 coefficients shown here 

1 

1 

r 13 s r l4 3 2 r 15 3 

i ne, * nj = g 

r 18 = r i 9 = r 20 = 0 

1 

r2 i s i 

1 

1 

r 22 3 r 23 * r 24 3 l 

r 25 3 0 r 26 3 | 

1 


r 27 =0 r 26 = - 

r 2 9 s r 3 o « 0 

1 

1 

r 31 3 1 r 32 3 l 

r 33 = r 3 j, * - r 3 5 


LOAD STATE VECTOR: 


r 1 « *1 

r 2 = X 2 

r 3 s X 3 

r« = X4 

r 5 3 x 5 

<2 

II 

X 

O' 


I 

input output I 

CALL "ACC" (r, . . . rg, k n , k 12 , k 13 ) I «• First functional evaluation 


For i s 1 to 3 ! r(i + 6) = r^ + ri 3 • aT • r(i ♦ 3 ) ♦ H5 • aT 2 • k^ 
r(i + 9) = r(i + 3) + r 2 6 ‘ AT • k 1A ; next i 
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Figure Cl.- Nystroo fourth order functional flow chart. 
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— 1 > 


Q 


input output 

CALL "ACC* (r^ . . . r-j 2 > k 2 , | , k 22 » k 23 ^ 


Second functional evaluation 


For i s I to 3 ; r(i + 6) s • AT • r(i + 3 ) + AT 2 (rj7 * k^ + ri8 * ) 

r(i + 9 ) = r(i + 3 ) + AT(r27 ’ k U + r 28 * k 2i> J next 1 


CALL "ACC" (r7 • . . r^ 2 i ^31 » k 32 » k 33 ^ 


Third functional evaluation 


For i : 1 to 3 I r(i + 6) = r* + * AT • r(i + 3 ) + AT 2 (rj9 • k^ + r20 • + r 2 i * ^3^) 

r(i + 9 ) - r(i + 3 ) + AT(r29 * k 1i ♦ r 30 * k 2i ♦ r 3i ' k 3i> • next 1 

1 

1 

1 

1 

CALL "ACC 1 * (r ? • . . r 12 » *in > k i<2> *43) I ♦ Fourth functional evaluation 


For i s 1 to 3 ; X* = r± + AT • r(i + 3 ) ♦ A T 2 (r 2 2 * k^ + ^3 ' k 2 i r 2 i* * k3* + r25 • ty*) 
X ( i ♦ 3 ) * r( i ♦ 3 ) ♦ AT(r32 * k*u + ^3 * k 2 i + ^4 • k3^ + ^5 * k^) ; next i 


/ \ 

/ \ 

No / \ Yes ! 

/TcURR =\ ^ 1 PRINT: FINAL STATE VECTOR 

\ T f / ! AND STOP 

\ / i 

\ / 

\ / 


Figure Cl.- Concluded. 
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APPENDIX D 

INTEGRATOR COEFFICIENTS 
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TABLE D1 . - INTEGRATOR COEFFICIENTS 



RK3 

RKL3 

NXD3 

NLXD3 

RK4 

r l3 

1/3 

6 -Vfa 
10 

1/2 

0.6 - \Tob 

1/2 

r 1 4 

2/3 

6 +Vb 
10 

1 

.6 + yfoe 

1/2 

r 15 

1/3 

6 -Vo" 

1/8 

.21 - .6^00 

1 

H6 

0 

-(54 + 19 V6)/250 

0 

( .15 + 4 V-06J/25 

1/2 

r 17 

2/3 

(102 ♦ 22 V6)/125 

1/2 

(5.1 + 11 ^0b)/25 

0 

r l8 

1/4 

1/9 

1/6 

1/9 

1/2 

r 1 9 

0 

(16 + V6)/36 

1/3 

(7 + 20 yp)b)/3(> 

0 

r 20 

3/4 

(16 -Vo)/36 

0 

(7 - 20 V^06)/36 

0 

r 21 



1/2 

(.6 - VTo6) 

1 

r 22 



-1 

-(5.4 ♦ 19 Vo6)/25 

1/6 

r 23 



2 

(20.4 ♦ 44 V^06)/25 

1/3 

r 24 



1/6 

1/9 

1/3 

r 25 


o c 
« 2 

ro c 
o ^ 

i; 
> - 

2/3 

(8*5 ^oT)/l8 

1/6 

r 26 

r 27 

r 28 

r 29 


1/6 

(8-5 y{^ot)/y& 


r 30 

r 3l 


1: 




r 32 







RKG4 

RKL14 

RKL24 

NXD4 

NLXD4 

1/2 

0.15 

0 

1 

(T> 

1/2 

(5 -V51/10 

1/2 

.192 

(5 + V^)/10 

1/2 

(5 ♦ Vi ) /10 

1 

1 

1 

1 

1 

1/2 

.15 

(5 -V5)/10 

1/8 

(3 - Vs) /20 

(-i *y[2)/2 

.1536 

-(5 ♦ 3 Vii /20 

1/8 

0 

(2 -V2)/ 2 

.0384 

(3 ♦ Vi>/4 

0 

(3 ♦ Vi) /20 

0 

6.74526571119 

(-1 ♦ 5 Vi)/4 

0 

(-1 ♦ Vi) /4 

-V2/2 

-38.7783195429 

-(5 ♦ 3 Vi)/4 

0 

0 

(2 *yfz )/2 

33.0330538317 

(5 -V5)/2 

1/2 

(3 - 5)/4 

1/6 

1 .41435185185 

1/12 

1/6 

1/12 

(2 -V2)/6 

-9.58605664488 

5/12 

1/6 

(5 ♦Vi)/24 

(2 +V 2 )/ 6 

8.95271818848 

5/12 

1/6 

(5 -Vi)/2 4 

1/6 

.218986604542 

1/12 

0 

0 




1/2 

(5 - Vii /10 




0 

-(5 ♦ 3Vi)/20 




1/2 

(3 *Vi)/4 




0 

(-1 ♦ sVi)/4 




0 

(-5 ♦ 3 Vi) /l4 




1 

(5 -Vi)/2 




1/6 

1/12 


r 
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TABLE DU- Concluded 

RK3 

RKL3 

NXD3 

NLXD3 

RK4 

RK G4 

RKi-14 

RKL24 

NXD4 

NLXDU 

r 33 








1/3 

5/12 

rjn 








1/3 

5/12 

r 35 








1/6 

1/12 


Note : 


In this study, only J2 is considered (except for reference 5U Therefore, for third order integrators, 
are not used. Similarly, for fourth order integrators, ri 3 , ri 4 , and r^ constants are not used. 


ri} and rj 4 


constants 
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APPENDIX E 

ACCELERATION FUNCTION SUBROUTINE 


E-1 
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21 

O 


! 


i 

CJ 


o 

>51 

G 

G 

£ 


o 

:) 

o 

> 

r 

> 

o 

P5 


*< 00 


C 


START 


J 


LOAD CONSTANTS: 

= 3.98601D5; J 2 = 1.08265D-3; ft B = 6.37122D3 
P DBN 1 8. 3631D-1; C Ay = 0.67385D-2 
< 

— i Intermediate Calculation 

r MAG = ^n 2 ♦ *' 2 2 + r 3 2 " ; 

1 

R 5IN = < R 7IN = 

R MAG 

X LAM * 1.5 • U B • J 2 * R B 2 
X LTB - X LAM ( ~5 * r 3 2 * r 7IN 


COMPUTE 

CENTRAL FORCE ACCELERATION: 


8CF = - 

V E 


r mag 



fbr i 

S 1 to 3 ; AcFi = b CF * r i > 

next i ! 


Vmag = v/v, 2 + V 2 2 ♦ V3 2 

R 7 

r mag 

+ r 5IN ) 
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Figure El.- Acceleration function subroutine 
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m 

i 

-p* 



COMPUTE J2 (OBLATENESS) ACCELERATION: 

for i : 1 to 2 ; Aj 2 i = ~ X LTB * r i » next i 

*J23 = -X LTB " r 3 _ X LAM ‘ 2 ‘ r 3 * R 5IN 



COMPUTE DRAG ACCELERATION: 


«B 

RE s 



Z = (R^q - Rg)/0. 30*48 I 

I 

I 

Z FUNCT = ® ( -24 .2 — 0* 00289 • Z ♦ 2605/Z) ! 

i 

DRAG = -0,5 • ^DBN * Z FUNCT * V MAG • 

i 

for i = 1 to 3 ; Ajji = DRAG ■ V A ; next i ! 

- - - 4 

I 

I 

I 

x 

TOTAL ACCELERATION: 

For i = 1 to 3 J ^TOTALi = Aq? ^ « Aj 2i «■ A^ ; 



Note 1: Stat^ vector units. 

A 

r 1 » r 2 > r 3 = position In km 
A 

V 1f V 2# V 3 = velocity in km/sec 


NOTE 2: r^ = position computation along 

Earth’s North Pole. 


next i 
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El.- Concluded. 



